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Abstract 


A  two-phase  study  was  conducted  to  determine  the  efficacy  of  MIL-STD-662E  in  balhstic 
protection  limit  determination  for  advanced  lightweight  armor  technologies.  The  &st-phase 
results  demonstrated  an  unacceptable  variability  in  limit  estimates  across  sites.  It  was  concluded 
that  more  stringent  controls  in  implementing  MIL-STD-662E,  particularly  as  to  test  velocity 
selection,  would  be  required  to  achieve  adequate  reproducibility.  An  unusual  pattern  in  the  data 
from  the  first  phase  suggested  that  there  was  a  physical  reason  beyond  operational  concerns  for 
the  unacceptable  variability.  In  the  second  phase,  this  new  reason  was  investigated  in  an 
interlaboratory  test  involving  six  sites.  It  was  determined  that  a  physical  phenomenon,  termed 
shatter  gap,  was  the  principal  cause  of  the  irreproducible  results  from  the  first  phase.  For  this 
setting,  modifications  of  MIL-STD-662E  to  support  reproducible  results  were  suggested  and 
successfully  tested. 
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1.  Introduction 


The  ballistic  performance  of  lightweight  armor  materials  has  been  evaluated  based  on  the  V50 
ballistic  limit,  which  is  defined  as  the  projectile  striking  velocity  at  which  complete  penetration  and 
partial  penetration  of  the  armor  materials  are  equally  likely  events.  There  are  several  types  of  Vjo 
ballistic  limits,  which  depending  on  the  definition  of  the  complete  (or  partial)  penetration,  included 
the  Army,  Navy,  and  protection  criterion  [1].  The  protection  ballistic  limit  (PEL)  Vjq  is  the  most 
widely  accepted  criterion  for  assessing  the  performance  of  lightweight  armor  materials  and  serves 
as  the  focus  of  this  study.  The  outcome  of  the  PEL  V5(,  baUistic  test  is  determined  by  the  final 
condition  of  a  wimess  plate  placed  behind  the  armor  panel.  If  the  wimess  plate  is  perforated  by  the 
projectile  or  spall  from  the  armor  panel,  the  outcome  is  termed  a  complete  penetration.  If  no 
perforation  is  observed  throu^  the  witness  plate,  the  result  is  termed  a  partial  penetration.  If  the 
armor  panel  is  perforated,  but  the  witness  plate  remains  intact,  the  result  is  still  defined  as  a  partial 
penetration.  A  schematic  definition  of  the  partial  and  complete  penetration  is  shown  in  Figure  1. 


Partial 

Penetration 

Armor  Witness 


witness  plate  is  intact 


Complete 

Penetration 

Armor  Witness 


witness  plate  is  not  intact 


Figure  1.  A  Schematic  Illustration  of  Partial  and  Complete  Penetrations. 


The  outcome  of  a  given  ballistic  event  (partial  or  complete  penetration)  at  a  given  velocity  V  is 
not  always  identical,  due  to  the  statistical  nature  of  the  various  material  failure  processes  involved 
in  the  ballistic  event.  Therefore,  it  is  best  described  in  terms  of  the  probability  of  complete 
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penetration,  Pc(V).  There  is  no  a  priori  derivation  of  Pc(V)  from  first  principles.  It  is  generally 
accepted  that  Pc(V)  is  assumed  to  follow  some  unspecified  statistical  model.  In  general,  Pc(V)  -  0 
at  lower  velocities,  Pc(V)  —  0.5  at  V50,  where  the  partial  and  complete  penetration  are  equally  likely, 
and  Pc(V)  -  1  at  higher  velocities.  Near  V50,  some  of  the  partial  penetrations  have  velocities  higher 
than  some  of  the  complete  penetrations.  The  zone  of  mixed  results  (ZMR)  is  the  velocity  range  over 
which  such  overltq)  occurs.  Ideally,  Pc(V)  can  be  determined  expaimentally  by  conducting  multiple 
ballistic  tests  at  various  velocities.  However,  this  approach  is  not  practical  for  obtaining  the  V50 
ballistic  limit  of  armor  materials  on  a  regular  basis,  due  to  the  high  cost  associated  with  such  large 
number  of  ballistic  tests.  Therefore,  a  more  efficient  and  standardized  test,  namely  MIL-STD-662E 
[2],  using  the  up-and-down  velocity-selection  rule  to  sample  the  Pc(V)  with  fewer  numbers  of 
ballistic  tests,  was  developed,  and  it  has  been  widely  used  throughout  the  ballistic  testing 
community.  The  up-and-down  method  is  based  on  a  bisection  algorithm.  If  die  nth  shot  at  V„  yields 
a  complete  penetration,  the  velocity  of  the  (n  +  l)th  shot  should  be  at  V„^.i  “  -  AV  in  order  to 

obtain  a  partial  penetration.  Otherwise,  if  the  nth  shot  yields  a  partial  penetration,  V^j  should  be  at 
V„ + AV  in  order  to  obtain  a  complete  penetration.  The  velocity  step,  AV,  is  typically  100  ft/s.  After 
the  occurrence  of  the  first  reversal,  that  is  a  complete  penetration  followed  by  a  partial  penetration 
or  a  partial  penetration  followed  by  a  complete  penetration,  AV  -  50  ft/s  and  AV  «  25  ft/s  for 
additional  reversals.  The  up-and-down  method  assumes  that  Pc(V)  is  monotonic;  (i.e.,  Pc(V)  is 
always  increasing  with  increasing  striking  velocity  in  the  velocity  range  of  interest). 

MIL^STD-662E  Ballistic  Test  for  Armor  [2]  defines  the  ballistic  limit,  protection  criteria,  Vjq 
BL(P),  as  follows. 

The  Vjo  BL(P)  may  be  defined  as  the  average  of  an  equal  number  of  highest  partial 
penetration  velocities  and  the  lowest  complete  penetration  velocities  that  occur 
within  a  specified  velocity  spread.  The  normal  up-and-down  firing  procedure  is  used. 

A  0.20-in  (0.051  mm)-thick  2024  T3  sheet  of  aluminum  is  placed  6  ±0.5  in  (152 
±12.7  mm)  behind  and  parallel  to  the  target  to  wimess  complete  penetrations. 
Normally,  at  least  two  partial  and  two  complete  penetration  velocities  are  used  to 
complete  the  BL(P).  Four-,  six-,  and  10-round  ballistic  limits  are  frequently  used. 
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The  mayimiim  allowable  velocity  span  is  dependent  on  the  armor  material  and  test 
conditions.  Maximum  velocity  spans  of  60,  90, 100,  and  125  ft/s  (18,  27,  30,  and 
38  m/s)  are  frequently  used. 

With  die  advancement  of  armor  material  technology,  especially  with  the  combinations  of  ceramic 
front  tile  with  composite  backing  plate,  there  have  been  many  instances  where  V50  data  obtained 
from  different  laboratories  differ  from  one  another.  Furthermore,  in  balhstic  tests  conducted  by  the 
same  laboratory,  there  were  instances  where  the  velocity  spreads  between  the  highest  partial 
penetration  and  the  lowest  complete  penetration  exceeded  the  allowable  limits.  In  order  to  validate 
the  cuiTMit  MI1^STD-662E  V50  test  procedure,  an  interlaboratory  ballistic  test  program  was  initiated 
by  the  U.S.  Army  Research  Laboratory  (ARL),  Aberdeen  Proving  Ground  (APG),  MD,  in  1993. 

The  objectives  of  this  project  were  to  (1)  evaluate  the  reproducibility  of  the  V50  data  using  the 
current  ME^STD-662E  by  conducting  mterlaboratory  ballistic  tests  on  benchmark  armor  materials 
and  (2)  propose  modifications  to  MIL-STD-662E.  This  project  was  sponsored  by  the  Joint  Techmcal 
Coordinating  Group  on  Aircraft  Survivability  (JTCG/AS).  The  JTCG/AS  project  identification 
number  is  V-3-A1 ,  Aircraft  Armor  Materials  Testing  and  Data  Reporting  Standards  Project.  This 
project  was  divided  into  two  phases  as  follows. 

Phase  I  (FY93-FY95):  Nine  laboratories,  including  four  government  laboratories  (ARUMD, 
Aviation  Applied  Technology  Directorate  [AATD],  U.S.  Air  Force  Wright  Laboratory  [AFWL],  and 
Naval  Surface  Warfare  Center  [NSWC]),  and  five  industrial  laboratories,  (H.  P.  White,  Ceradyne, 
Simula,  Southwest  Research  Institute  [SRI],  and  University  of  Dayton  Research  Institute  [UDRI]), 
participated  in  this  project.  These  laboratories  conducted  Vjq  ballistic  tests  in  accordance  with 
MIL^STD-662E  [2]  on  identical  armor  materials  (material  no.  1:  VAR  4340  Steel,  material  no.  2: 
AD90  AljOj/Kevlar)  using  identical  0.50  AP  M2  bullets  from  the  same  lot.  Results  from  Phase  I 
[3]  demonstrated  the  weakness  of  the  MIL-STD-662E  procedure  for  obtaining  consistent  Vjo  data 
from  different  laboratories.  Statistical  analysis  of  Phase  I  data  suggested  that  Pc(V)  for  both  4340 
steel  and  AljOj/Kevlar  armor  plates  were  nonmonotonic.  This  nonmonotonic  behavior  implies  that 
each  laboratory  might  have  sampled  different  parts  of  this  Pc(V)  by  using  slightly  different  starting 


3 


velocities,  thus  obtaining  significantly  different  V5Q  data.  Other  variations  in  experimental 
techniques  from  each  laboratory,  such  as  different  fixtures  for  mechanically  constraining  the  armor 
plates  during  impact,  would  increase  the  data  scatter  in  the  interlaboratory  test  even  further. 

Phase  n  (FY95-FY96):  A  new  V50  test  procedure  was  developed  during  this  phase.  This 
proposed  V50  test  procedure  can  sample  both  the  monotonic  and  the  nonmonotCMiic  Pc(V)  and  was 
successfully  verified  in  Phase  n  interlaboratory  V50  ballistic  tests  that  involved  six  laboratories, 
including  three  government  laboratories  (ARL/MD,  AATD,  and  AFWL)  and  three  industrial 
laboratories  (Ceradyne,  Simula,  and  SRI).  The  Phase  n  results  are  the  main  subject  of  this  report. 

In  the  following  sections,  we  briefly  summarize  the  Phase  I  results  and  identify  the  sources  of 
discrepancy  in  Phase  I  data,  describe  efforts  in  confirming  the  nonmonotonic  behavior  of  Pc(V)  in 
AljOj/Kevlar  armor  plates,  and  describe  the  development  of  the  new  velocity-selection  rule  for  the 
V50  test  by  performing  statistical  model  simulations.  This  new  velocity-selection  rale  was  field- 
tested  by  the  six  laboratories  in  the  Phase  n  interlaboratory  V50  ballistic  tests,  yielding  consistent  Vjq 
data  from  the  participants.  We  also  discuss  die  physical  processes  responsible  for  the  nonmonotonic 
behavior  and  evaluate  the  experimental  techniques  that  contributed  to  the  V50  data  scatter.  Finally, 
we  provide  a  proposed  V50  test  procedure. 

2.  Review  of  Phase  I  Results 

Two  types  of  armor  materials/systems  were  tested  by  nine  laboratories  in  Phase  I.  The  first  type 
is  the  0.175-in-thick  VAR  4340  steel  panels  with  a  Rockwell  C  hardness  of  48-52.  The  second  type 
is  the  ceramic/composite  panels  with  AD90  AI2O3  ceramic  tile  (5  in  x  5  in  x  0.535  in)  backed  by 
Kevlar-reinforced  plastic  (KRP)  (15  in  x  15  in  x  18  ply).  Each  laboratory  was  given  3  VAR  4340 
steel  panels  capable  of  sustaining  8  or  9  shots  each,  16  AljO/Kevlar  panels  capable  of  sustaining 
a  single  shot,  and  enough  0.50-cal.  AP  M2  bullets  that  came  from  the  same  lot  The  starting  velocity 
of  1,600  ft/s  was  selected  for  both  the  steel  panels  and  the  ceramic/composite  panels.  The  test 
engineer  of  each  laboratory  selected  all  subsequent  velocities  (by  selecting  the  powder  charges) 
according  to  the  velocity-selection  rale  specified  by  MIL-STD-662E.  All  laboratories  used  variants 
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of  the  up-and-down  method  to  select  the  projectile  velocities.  Most  of  the  participants  focused  their 
attention  in  the  velocity  range  from  1,400  to  2,000  ft/s  for  both  targets.  After  the  V50  PBL  had  been 
satisfactorily  determined,  each  laboratory  used  the  remaining  shots  (usually  two-four)  to  explore 
velocities  outside  the  1,400-2,000-ft/s  range. 

2.1  Ballistic  Data  and  Penetration  ProbabOity.  Figure  2  summarizes  the  ballistic  data  of 
VAR  4340  steel  armor  plates  from  all  nine  laboratories.  It  is  clear  that  the  velocity  range  where  the 
partial  and  complete  penetrations  overlap,  ZMR  varies  from  laboratory  to  laboratory.  For  example, 
in  Lab  6,  the  ZMR  was  as  large  as  600  ft/s,  while  in  Lab  4,  it  was  only  about  150  ft/s.  Therefore,  the 
reported  V50  data  varied  from  laboratory  to  laboratory  with  little  consistency.  Figure  3  shows  the 
ballistic  data  of  AljOa/Kevlar  armor  plates  from  all  nine  laboratories.  Again,  the  ZMR  varies  from 
laboratory  to  laboratory  and  appears  to  be  as  large  as  800  ft/s  in  Lab  9.  By  grouping  the  combined 

from  all  ninft  laboratories  into  various  velocity  bins,  the  penetration  probability  for  a  given 
velocity  bin,  Pc(V),  can  be  estimated  by  dividing  the  number  of  complete  penetrations  by  the  total 
number  of  shots  in  this  velocity  bin.  Figure  4  shows  the  Pc(V)  obtained  from  the  combined  data  of 
VAR  4340  steel  armor  plates,  and  Figure  5  shows  the  Pc(V)  obtained  from  the  combined  data  of 
Al203/Kevlar  armor  plates.  The  solid  squares  in  Figures  4  and  5  are  observed  penetration  probability 
at  various  velocities.  The  solid  curves  are  fitted  to  the  data  based  on  a  proposed  statistical  model, 
which  will  be  discussed  in  detail  later.  The  observed  penetration  probability  data  in  Figures  4  and 
5  strongly  suggest  that  the  Pc(V)  in  fliese  armor  materials  systems  is  nonmonotonic.  In  other  words, 
the  probability  of  complete  penetration  does  not  always  increase  with  increasing  velocity.  In  certain 
velocity  ranges,  Pc(V)  actually  decreases  with  increasing  velocity. 

2.2  Sources  of  Discrepancy.  There  are  two  kinds  of  discrepancy  in  tiiis  interlaboratory  ballistic 
test,  namely  the  extrinsic  discrepancy  vs.  the  intrinsic  discrepancy.  The  extrinsic  discrepancy  could 
come  from  the  difference  in  V50  ballistic  test  practice  in  each  laboratory  (such  as  difference  in  target 
fixture,  velocity  measurement  and  control,  bullet  pitch  and  yaw,  and  velocity  selection).  The 
intrinsic  discrepancy,  on  the  other  hand,  could  come  from  the  variation  of  material  behaviors  in 
bullets  and  armor  panels  as  a  function  of  bullet-striking  velocity.  This  intrinsic  discrepancy,  if  exists, 
should  be  observable  even  within  each  laboratory’s  own  data.  Careful  examination  of  Figures  2  and 
3  indicates  that  even  within  a  single  laboratory’s  data,  the  indication  of  a  nonmonotonic  Pc(V) 
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Laboratories. 


Lab 


Figure  3.  Summary  of  Phase  I  Ballistic  Data  of  AljOj  /KeVlar  Armor  Plates  From  All  Nine 
Laboratories. 
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Impact  Velocity  (ft/s) 


Figure  4.  Penetration  Probability  of  VAR  4340  Steel  Armor  Plates  Obtained  From  Combined 
Phase  I  Data  of  All  Nine  Laboratories. 


Impact  Velocity  (ft/s) 


Figure  5.  Penetration  Probability  of  AljOj  /Kevlar  Armor  Plates  Obtained  From  Combined 
Phase  I  Data  of  Ail  Nine  Laboratories. 
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distribution  is  evident.  Assuming  the  existence  of  a  nonmonotonic  Pc(V)  distribution,  it  is 
postulated  that  the  current  velocity-selection  rule  could  make  each  laboratory  sample  different  parts 
of  a  nonmonotonic  Pc(V)  distribution.  This  intrinsic  discrepancy,  if  not  corrected,  could  lead  to 
drastically  different  Vjq  results,  overwhelming  the  effect  of  the  extrinsic  discrepancy. 

Therefore,  the  strategies  of  the  Phase  n  study  were  (1)  to  confirm  the  existence  of  a 
nonmonotonic  penetration  probability  distribution,  Pc(V)  in  Al203/Kevlar  armor  plates  by 
conducting  systematic  ballistic  tests  at  various  velocities  by  a  single  laboratory,  thus  eliminating  the 
possible  effects  of  extrinsic  discrepancy  due  to  interlaboratory  inconsistencies  and  bypassing  the 
effect  of  the  current  velocity-selection  rule;  (2)  to  develop  a  new  velocity-selection  rule,  suitable  for 
both  the  nomnonotonic  and  the  monotonic  Pc(V)  distributions,  fi’om  statistical  model  simulations; 
and  (3)  to  confirm  the  benefit  of  the  new  velocity-selection  rule  by  conducting  the  Phase  n 
interlaboratory  V50  test. 

3.  Phase  n  Results 

In  the  following  sections,  we  present  the  experimental  and  analytical  details  of  the  confirmation 
of  the  nonmonotonic  behavior  of  Pc(V),  the  direct  observation  of  bullet  shattering  as  a  function  of 
velocity,  physical  evidence  of  the  shatter-gap  behavior,  the  development  of  statistical  models  to 
describe  this  nonmonotonic  Pc(V),  the  development  of  a  new  velocity-selection  mle  by  Monte-Carlo 
simulation  of  the  Vjo  test  procedure,  and  the  results  of  the  Phase  n  interlaboratory  V50  test  (involving 
six  laboratories)  using  this  new  velocity-selection  rule  for  nonmonotonic  Pc(V). 

3.1  Confirmation  of  Nonmonotonic  Penetration  Probability  Distribution.  A  new  batch  of 
144  AljOj/Kevlar  armor  plates,  identical  to  the  Phase  I  specifications,  was  procured,  hi  collaboration 
with  the  U.S.  Army  Aberdeen  Test  Center  (ATC),  ARL/MD,  conducted  69  ballistic  shots  to  identify 
the  Pc(V).  The  remaining  75  armor  plates  were  used  by  five  additional  laboratories  to  confirm  the 
new  Vjo  test  procedure.  Using  Phase  I  results  as  the  guide,  the  velocity  window  of  interest  was 
identified  to  be  from  1,400  to  2,800  ft/s.  The  velocity  window  was  divided  into  11  bins.  The 
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number  of  shots  per  bin  was  selected  depending  on  the  velocity.  For  velocity  bins  at  which 
nonmonotonic  behavior  was  anticipated,  a  greater  number  of  shots  were  allocated  to  increase  the 
confidence  of  die  observations.  Flash  x-ray  radiographs  and  yaw  cards  were  used  to  confirm  that 
the  0.50-cal.  AP  M2  bullet  had  acceptable  pitch  and  yaw  (less  than  3°). 

Figure  6  shows  the  penetration  probability  derived  from  the  69  shots.  The  solid  squares  in 
Figure  6  are  observed  penetration  probabilities  at  various  velocities.  The  upper  and  lower  80% 
confidence  intervals  associated  with  each  observed  probability  are  also  shown  as  hollowed  triangles 
and  hollowed  circles,  respectively.  The  solid  curve  is  based  on  a  data-fitted  statistical  model  that 
will  be  discussed  later. 


impact  Velocity  (ft/s) 


Figure  6.  Penetration  Probability  of  AljOj  /Kevlar  Armor  Plates  Obtained  From  69  Shots  in 
Phase  n. 

Table  1  lists  the  results  of  the  69  shots  along  with  the  calculated  Pc(V)  and  80%  and  90% 
confidence  intervals.  The  Pc(V)  of  a  velocity  bin  at  Vmean  is  simply  the  number  of  complete 
penetrations  (C),  divided  by  the  total  number  of  shots  in  this  velocity  bin. 
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Table  1.  Results  of  69  Ballistic  Shots  on  Phase  H  AljOj  /Kevlar  Armor  Plates 


Note:  P  -  Partial 
C  -  Complete 
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Table  1.  Results  of  69  Ballistic  Shots  on  Phase  H  Armor  Plates  (continued) 


P/C 

Pc(V) 

Lower  80% 

Upper  80% 

Lower  90% 

Upper  90%  I 

B^l 

mmM 

1,897 

C 

1,912 

4/4 

0.562 

1 

0.473 

1 

C 

c 

1,929 

c 

1,954 

p 

2,004 

3/6 

0.201 

0.799 

0.153 

0.847 

1,966 

c 

2,005 

p 

2,019 

p 

2,023 

c 

2,057 

c 

2,197 

p 

2,235.4 

0/5 

0 

0.369 

0 

0.451 

2,224 

p 

2,240 

p 

2,253 

p 

2,263 

p 

2,383 

p 

2,426.5 

1/8 

0.013 

0.406 

0.006 

0.471 

2,398 

p 

2,403 

p 

2,425 

p 

2,426 

p 

2,437 

p 

2,454 

2,486 

2,597 

c 

2,618.14 

6/7 

0.547 

0.985 

0.479 

0.993 

2,605 

c 

2,618 

p 

2,621 

c 

2,626 

c 

2,626 

c 

2,634 

c 

2,747 

c 

2,807.2 

5/5 

0.631 

1 

0.549 

1 

2,816 

c 

2,823 

c 

2,824 

c 

2,826 

c 

Note:  P  -  Partial 
C  -  Complete 
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3^  Observation  of  Shatter  Gap.  Rash  x-ray  was  also  used  to  record  and  measure  residual 
bullet  length  from  complete  penetrations.  For  partial  penetrations,  the  residual  bullet  was  collected 
and  the  residual  bullet  length  measured,  whenever  possible.  Based  on  the  residual  bullet  length 
measurement  from  both  the  partial  and  the  complete  penetrations,  the  bullet-shattraing  phenomenon, 
the  shatter  gap,  was  directly  confirmed  by  this  test.  Figure  7  shows  the  shattering  of  0.50-cal.  AP 
M2  buUet  after  impacting  AljOj/Kevlar  armor  plates  as  a  function  of  impact  velocity.  The  dotted 
curve  is  the  same  data-fitted  model  Pc(V)  as  shown  in  Figure  6.  Residual  core  lengths  obtained  from 
flie  partial  and  the  complete  penetrations  are  plotted  separately  as  hollow  squares  and  solid  triangles, 
respectively. 

2 
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0.2 
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0  500  1000  1500  2000  2500  3000 
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Figure  7.  Shattering  of  0.50-cal.  AP  M2  BuUets  as  a  Function  of  Velocity  After  Impacting 
AI2O3  /Kevlar  Armor  Plates. 

At  velocities  below  1,700  ft/s  where  partial  penetrations  were  most  probable,  the  residual  core 
length  reduced  only  slightly  from  tiie  original  1 .8  in  by  breaking  the  nose  tip.  At  velocities  between 
1,700  and  2,000  ft/s,  more  bullet  shattering  (cote  length  reduction)  was  observed  with  increasing 
velocity,  and  the  probability  of  complete  penetration  began  to  decrease  with  increasing  velocity. 
Between  2,000  and  2,400  ft/s,  where  no  complete  penetrations  were  observed,  we  did  not  recover 
any  residual  bullet  in  partial  penetrations.  At  this  velocity  range,  it  was  believed  that  the  bullet  was 
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so  severely  shattered  that  it  was  extremely  difficult  to  recover  the  residual  bullet  from  partial 
penetrations.  At  velocities  between  2,400  and  2,800  ft/s,  the  residual  core  length  showed  an 
increase,  probably  due  to  the  large  deformation  of  the  Kevlar  backing  plate  reducing  the  impact 

stress  on  the  bullet  The  probability  of  complete  penetration  increased  witii  increasing  velocity.  This 
observed  shatter-gap  behavior  is  believed  to  be  the  physical  processes  responsible  for  the  observed 
nonmonotonic  behavior  of  Pc(V). 

33  Statistical  Models  and  Simulation.  Adams  and  Zwissler  [4]  suggested  that  the 
nonmonotonic  Pc(V)  can  be  expressed  in  terms  of  a  combination  of  three  probability  distributions 
such  that 


Pc(V)  -  (1  -  P(Tm,  Ts,  V))  ♦  P(Im,  Is,  V)  +  P(Tm,  Ts,  V)  ♦  P(Sm,  Ss,  V),  (1) 


where 


P(Im,  Is,  V)  -  probability  of  complete  penetration  by  an  intact  bullet  at  V; 

P(Sm,  Ss,  V)  -  probability  of  complete  penetration  by  a  shattered  bullet  at  V; 

P(Tm,  Ts,  V)  -  probability  of  transition  between  intact  (0.0)  and  shatter  (1.0)  at  V; 

Im,  Sm,  Tm  -  mean  of  intact,  shatter,  and  transition,  respectively;  and 

Is,  Ss,  Ts  “  standard  deviation  of  intact,  shatter,  and  transition,  respectively. 

Equation  (1)  defines  Pc(V)  as  a  sum  of  two  probability  contributions.  The  first  one  is  from 
having  an  intact  bullet  at  V  that  completely  penetrates  the  target,  the  (1  -  P(Tm,  Ts,  V))  *  P(Im,  Is, 
V)  term.  The  second  one  is  from  having  a  shattered  bullet  at  V  that  completely  penetrates  the  target, 
the  P(Tm,  Ts,  V)  ♦  P(Sm,  Ss,  V)  term.  It  was  found  that  this  model  can  adequately  describe  the 
observed  nonmonotonic  Pc(V)  with  appropriately  chosen  parameters.  Table  2  shows  the  data-fitted 
parameters  used  to  describe  the  observed  Pc(V)  of  VAR  4340  Steel,  Phase  I  Al203/Kevlar,  and 
Phase  n  AljOj/Kevlar  armor  plates. 
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Table  2.  Data-Fitted  Model  Parameters  From  Observed  Pc(V) 


Armor  Plates 

B9 

Is 

(ft/s) 

— 

Ts 

Ss 

(ft/s) 

VAR  4340  Steel 

1,425 

125 

1,700 

125 

1,825 

125 

Phase  I  AljOs/Kevlar 

1,750 

125 

2,000 

125 

2,500 

125 

Phase  n  Al,0^evlar 

1,650 

100 

2,050 

100 

2,550 

100 

3^  New  Velocity  Selection  Rule.  The  current  program  offered  a  unique  opportunity  to 
re-examine  various  assumptions  of  the  MIL-STD-662E  V50  test  procedure.  The  MIL-STD-662E 
selection  of  initial  velocity  at  the  “expected  V50”  and  the  choice  of  reducing  the  step  size  after  an 
observed  reversal  were  based  on  the  assumption  of  a  monotonic  Pc(V),  such  that  the  subsequent 
velocities  would  converge  rapidly  toward  the  “actual  V50.”  The  MIL-STD-662E  method  of  V50 
calculation  by  averaging  equal  numbers  of  partial  and  complete  penetrations  around  the  ZMR  was 
also  based  on  the  assumption  that  Pc(V)  is  always  monotonic  and  Pc(V)  is  symmetrical  around  V5Q. 
Since  the  assumption  of  a  monotonic  Pc(V)  was  nullified  by  the  direct  measurement  of  a 
nonmonotonic  Pc(V)  in  Al203/Kevlar  systems,  modifications  to  various  aspects  of  the  Vjq  test 
procedure  were  studied  in  detail  with  Monte-Carlo  simulations. 

It  was  immediately  evident  that  in  a  system  with  nonmonotonic  Pc(V),  the  value  of  V50  was 
dictated  by  the  initial  velocity  selection.  Therefore,  the  only  sensible  alternative  was  to  have  two 
starting  velocities,  one  firom  the  low  end  of  the  velocity  window  of  interest  and  going  up  and  the 
other  from  the  high  end  and  coming  down.  In  othw  words,  one  should  perform  two  independent  V50 
tests  and  report  both  values  of  low- Vjq  and  higJi-Vjo.  If  the  two  values  agreed  within  the  uncertainty 
limits,  the  system  had  a  monotonic  Pc(V)  and  the  V50  -  (I0W-V50  +  high-Vjo)/2.  Otherwise,  the 
system  was  nonmonotonic,  and  one  should  report  both  values.  From  the  personnel  protection  point 
of  view,  the  low-Vj^  value  should  be  more  significant;  therefore,  a  greater  number  of  shots  should 
be  aUocated  toward  the  low-  V50  test  Compromise  between  statistical  accuracy  and  conservation  of 
resources  led  to  the  allocation  of  nine  shots  for  the  low  -V50  test  and  six  shots  for  the  high- V50  test 
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Important  to  both  design  and  estimation  is  the  selection  of  velocity  step  size.  In  the  design,  it 
governs  what  data  are  collected.  In  estimation,  it  impacts  the  V50  estimate.  One  danger  of  selecting 
a  step  siye  too  small  is  that  if  Pc(V)  is  truly  monotonic,  the  modest  sample  size  (nine  or  six)  might 
not  pomit  the  design  to  travel  quickly  enough  to  the  lone  Vjq.  When  the  step  size  is  too  large,  the 
sequential  design  “walks”  past  the  I0W-V50.  This  causes  the  estimates  for  V50  to  be  more  variable 
and  biased  toward  the  right. 

Although  by  no  means  rigorous,  we  would  suggest  discussions  along  the  following  lines  to  set 
tile  step  size,  we  suggest  that  the  tester  first  guesses  the  lowest  velocity,  Vi^^,  at  which  a  complete 
penetration  would  be  observed  and  then  the  highest  practicable  velocity,  Vhigi,,  based  on  physical 
considerations  of  the  gun,  bullet,  etc.,  at  which  a  partial  penetration  would  be  observed.  and 
Vi,i^  would  define  the  velocity  window  of  interest.  The  velocity  window  of  interest  for  the  current 
system  was  determined  from  Phase  I  data  to  be  between  -  1,400  and  -  2,800  ft/s.  We  also 
suggest  that  more  shots  be  taken  for  the  I0W-V50  than  for  high-VjQ.  In  this  project,  we  used  nine 
shots  for  tile  low-Vjo  -  9)  and  six  shots  for  the  high-Vjo  -  6).  In  order  to  make  sure  that 
both  V  tests  can  sample  a  monotonic  Pc(V)  situated  somewhere  in  the  middle  of  the  velocity 
window,  we  suggest  that  the  selected  step  size  would  provide  an  overlap,  AVo^^jjp,  in  the  two 
sampling  ranges.  Therefore,  the  step  size  AV  can  be  roughly  estimated  to  be 

AV  -  (V^gH  -  +  N^gh),  (2) 

In  the  current  test,  we  have  AV  -  (2,800  -  1,400  +  100)/(9  +  6)  -  100  ft/s. 

The  intended  velocity  for  the  first  shot  of  the  I0W-V50  test  was  set  at  Vj  -  +  AV,  and  the 

subsequent  eight  intended  velocities  were  set  at  Vj+i  -  Vj  +  AV,  if  the  previous  shot  showed  a  partial 
penetration  and  Vki-V,  -  AV  for  complete  penetrations.  Similarly,  the  intended  velocity  for  the 
first  shot  of  the  high- Vj^  test  was  set  at  Vj^  -  V^^  -  A  V,  and  the  subsequrait  five  intended  velocities 
were  set  at  Vj^,  -  Vj  -h  AV  for  partial  penetrations  and  Vj+j  -  Vj  -  AV  for  complete  penetrations. 
This  is  a  modified  up-and-down  method  without  reducing  the  step  size  after  an  observed  reversal. 
The  reason  for  using  the  intended  velocities  (instead  of  the  actual  velocities)  for  selecting  the  next 
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velocity  was  to  average  out  the  fluctuations  in  velocity  control  (the  gun  noise)  from  the  actual 
velocity  sequence.  In  simulations,  flie  gun  noise  was  assumed  to  have  a  maximum  range  of  ±25  ft/s, 
and  the  random  velocity  fluctuations  were  assumed  to  be  centered  around  each  of  the  intended 
velocities. 

The  effect  of  velocity  step  size  and  of  up-and-down  velocity  selection  on  the  simulated  V50 
distribution  was  studied  in  detail.  Figure  8  shows  the  effect  of  the  velocity  st^  size  (AV  “  50, 100, 
and  200  ft/s)  on  the  simulated  low-Vjo  distribution,  with  the  starting  velocity  Vj  -  1,400  +  50; 
1 ,400  + 100;  and  1,400  +  200  ft/s.  Figure  9  shows  the  same  for  the  simulated  high-Vjo  distribution, 
with  the  starting  velocity  Vjq  “  2,800  -  50;  2,800  -  100;  and  2,800  -  200  ft/s.  These  results 
indicated  that  the  choice  of  step  size  betwewi  50  and  (through)  200  ft/s  did  not  change  the  V50 
distribution  significantly. 


Velocity  (ft/s) 


Figure  8.  Effect  of  Velodty  Step  Size  on  the  Simulated  Low-Vjo  Distribution  in  Armor  Plates 
With  a  Nonmonotonic  Pc(V). 

Figure  10  shows  that  the  proposed  V50  test  procedure  is  also  applicable  to  armor  plates  with  an 
assumed  monotonic  Pc(V).  The  overlapping  low-Vjo  and  high-Vjo  distributions  (with  a  step  size  of 
100  ft/s)  suggested  that  the  actual  V50  -  (low- Vjq  +  high-Vjo)/2. 
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Figure  9.  Effect  of  Velocity  Step  Size  on  the  Simulated  High'Vso  Distribution  in  Armor  Plates 
With  a  Nonmonotonic  Pc(V)« 
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Figure  10.  Overlapping  of  Low  -Vjo  and  High-Vjo  Distributions  in  Armor  Plates  With  a 
Monotonic  Pc(V). 
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The  proposed  Vjo  test  procedure  is  summarized  in  Appendix  A.  This  procedure  was  sent  out  to 
the  participating  laboratories  for  the  Phase  II  interlaboratory  V50  test  program. 

3^  Modified  Dixon’s  Maximum  Likelihood  Estimate.  Since  there  is  no  guarantee  that  Pc(V) 
is  symmetrical  around  ZMR  in  a  nonmonotonic  Pc(V),  averaging  equal  numbers  of  partial  and 
complete  penetrations,  as  suggested  by  MIL-STD-662E,  would  not  always  give  the  best  estimate  of 
the  Vjo  value.  Instead,  we  have  developed  a  modified  Dixon’s  method  for  calculating  the  maximum 
likelihood  estimate  (MLE)  of  V50  (known  as  MDMLE). 

The  MLE  method  selects  a  V50  from  a  data  set  of  complete  and  partial  penetrations  with  their 
respective  velocities  that  would  make  V50  most  likely  to  have  occurred.  Estimation  of  I0W-V50  and 
high-V5o  from  two  separate  sets  of  data  was  accomplished  through  a  modification  to  the  usual  MLE 
approach.  Following  Dixon’s  [5]  approach,  the  standard  deviation  of  the  response  curve  was 
assumed  to  be  known  and  equal  to  the  fixed  step  size  taken  in  the  design.  Under  this  modification, 
the  conditions  for  estimate  existence  are  only  that  at  one  complete  penetration  and  one  partial 
penetration  are  observed  for  each  of  the  low-Vjo  and  high-VjQ  test  A  logistic  response  function  was 
selected  to  locally  describe  the  response  curve  about  each  Vjo-  The  Newton-Raphson  approximation 
method  was  used  to  cany  out  this  estimation  procedure.  A  FORTRAN  code  listing  of  program 
MDMLE  is  included  in  Appendix  B. 

3.6  Phase  n  Interlahoratory  Vg^  Test  hi  the  Phase  n  interlaboratory  V50  test  program,  each 
laboratory  followed  the  proposed  V50  test  procedure  (Appendix  A),  instead  of  the  MII^STD-662E 
procedure.  In  addition  to  identical  armor  plates  and  bullets,  the  same  identical  target  fixture  was 
used  by  all  six  laboratories.  The  reported  low-Vjo  and  high-Vjo  values  from  each  laboratory  were 
calculated  by  ARL/MD  using  the  identical  MDMLE  program.  Table  3  lists  the  results  of  the  Phase  n 
interlaboratory  V50  test  from  six  laboratories.  Figure  1 1  is  the  graphical  representation  of  Table  3. 
The  measured  low-  V50  and  high-  V50  values  from  each  laboratory  are  shown  as  solid  squares  with  an 
‘x.’  The  dotted  lines  represent  the  Iow-Vjq  and  high-Vjo  values  obtained  from  the  model  V50  values 
fitted  to  the  observed  Pc(V).  The  consistency  among  the  interlaboratory  V50  data  is  excellent,  far 
superior  than  that  obtained  in  Phase  I  using  MIL-STD-662E.  Figure  12  displays  the  observed  V50 
values  with  simulated  distributions  for  these  values  under  the  model  [1]. 
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Table  3.  Results  of  Vj#  Ballistic  Test  on  Phase  n  AljOj  /Kevlar  Armor  Plates  From  all  Six 

Laboratories 


Lab 

Velocity 

P/C 

L0W-V50 

Lab 

Velocity 

P/C 

High-Vjo 

1 

P 

1,609 

1 

2,747 

C 

2,475 

1 

C 

1 

2,597 

C 

1 

1,518 

1 

2,454 

P 

1 

1,597 

1 

2,621 

c 

1 

1,679 

1 

2,486 

c 

1 

1,609 

1 

2,403 

P 

1 

1,705 

1 

1,616 

1 

1,492 

2 

m 

P 

1,727 

2 

2,692 

c 

2,524 

2 

B& 

P 

2 

2,504 

p 

2 

Isil 

C 

2 

2,596 

c 

2 

■i9 

C 

2 

2,492 

c 

2 

P 

2 

2,406 

p 

2 

P 

2 

2,500 

p 

2 

BiS 

P 

2 

C 

2 

KB 

C 

3 

P 

1,647 

3 

2,712 

c 

2,543 

3 

C 

3 

2,600 

c 

3 

P 

3 

2,511 

p 

3 

1,604 

P 

3 

2,608 

c 

3 

1,730 

C 

3 

2,532 

p 

3 

1,604 

C 

3 

2,601 

c 

3 

1,474 

P 

3 

1,609 

P 

3 

1,703 

P 

4 

P 

1,632 

4 

2,732 

c 

2,570 

4 

C 

4 

2,623 

p 

4 

P 

4 

2,680 

c 

4 

1,646 

C 

4 

2,597 

c 

4 

1,447 

P 

4 

2,495 

p 

4 

1,648 

P 

4 

2,597 

c 

4 

1,738 

C 

4 

1,637 

P 

4 

1,656 

C 

Note:  P  -  Partial 
C  -  Complete 
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Table  3.  Results  of  V50  Ballistic  Test  on  Phase  n  AljOj  /Kevlar  Armor  Plates  From  All  Six 
Laboratories  (continued) 


Note:  P  -  Partial 
C  -  Conq)lete 

Results  of  this  test  confirmed  the  usefulness  of  the  proposed  Vjq  test  procedure  and  provided  the 
valuable  operational  experience  for  future  implementation  of  this  procedure  at  large. 

4.  Discussion 


The  current  program  had  the  benefit  of  a  well-characterized  Pc(V),  which  directly  confirmed  the 
shatter-gap  behavior  postulated  by  many  previous  investigators.  Armed  with  a  statistical  model  of 
this  nonmonotonic  Pc(V)  and  a  Monte-Carlo  simulation  program,  we  can  verify  various  strategies 
for  V50  test  and  quantify  the  build-in  uncertainties  in  V50  test. 


Lab 


Figure  11.  Summary  of  Phase  IT  Vjq  Ballistic  Test  on  AI3O3  /Kevlar  Armor  Plates  From  All  Six 
Laboratories. 


Figure  12.  Measured  Low-Vjo  and  High-Vjo  From  Phase  n  Interlaboratory  Test  on 
AI2O3  /Kevlaf  Armor  Plates  as  Compared  to  Simulations. 
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4.1  Shatter  Gap.  The  exact  physical-loading  conditions  and  mechanisms  for  initiating 
extensive  bullet  shattering  are  not  known  at  this  point  More  instramented  ballistic  tests  are  needed. 
It  was  noted  that  as  the  bullet  shattering  progressed  with  increasing  velocity,  the  Kevlar  backing  plate 
began  to  deform  significantly,  forming  a  big  bulge  without  perforation  for  partial  penetrations.  The 
large  deformation  of  the  backing  plate  seemed  to  decrease  the  loading  on  the  bullet  and  reduce  the 
bullet  shattering.  At  even  higher  velocities,  more  complete  penetrations  began  to  appear  and  the 
severely  deformed  backing  plate  began  to  form  a  hole  by  the  residual  bullet 

The  existence  of  a  nonmonotonic  Pc(V)  in  armor  systems  has  profound  implications  in  regard 
to  the  armor  system  acceptance  tests.  Lack  of  complete  penetrations  at  muzzle  velocities  does  not 
guarantee  that  the  system  will  not  be  penetrated  at  lower  velocities.  It  is  hoped  that  the  proposed  V50 
test  procedure  will  serve  as  the  basis  for  a  industry-adopted  test  standard. 

4.2  Vjo  Data  Scatter.  Due  to  the  small  number  of  shots  allocated  and  the  build-in  ±25  ft/s  gun 
noise,  the  simulated  I0W-V50  distribution  (nine  shots,  AV  -  100  ft/s)  was  centered  at  1,650  ft/s  with 
a  spread  width  of  ~120  ft/s  (measured  between  the  10th  and  the  90th  quantile)  while  the  simulated 
high- Vjo  distribution  (six  shots,  AV  -  100  ft/s)  was  centered  at  2,550  ft/s  and  had  a  spread  width  of 
~140  ft/s.  These  spread  widths  represented  the  best  accuracy  in  V50  values  obtainable  under  ideal 
conditions.  In  reality,  the  velocity  measurement  apparatus  and  technique  of  various  laboratories 
were  seldom  calibrated  with  each  other.  Some  laboratories  applied  air-drag-correction  procedures 
to  the  measured  velocity,  while  others  did  not  Few  laboratories  used  flash  x-ray  to  monitor  bullet 
pitch  and  yaw.  Most  used  yaw  cards  instead  to  cut  down  cost  It  was  found  that  both  ends  of  the  gun 
barrel  must  be  rigidly  restrained  to  obtain  consistently  low  pitch  and  yaw.  In  the  current  program, 
all  six  laboratories  used  the  same  target  fixture  that  was  used  to  generate  the  Pc(V).  In  future 
implementations,  difference  in  target  fixture  geometry  and  attachment  technique  from  each 
laboratory  will  undoubtedly  contribute  to  additional  data  scatter.  Perfiaps  the  most  significant  source 
of  interlaboratory  data  scatter  came  from  the  large-gun  noise,  which  well  exceeded  die  ±25  ft/s  value 
used  in  simulations.  Laboratories  situated  in  dry-climate  regions  tend  to  get  lower  gun  noise  than 
those  located  in  humid-climate  regions. 
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5.  Conclusions 


Specifically,  we  conclude  the  following. 

(1)  The  cunent  MIL-STD-662E  V50  test  procedure  was  based  on  the  assumption  that  the  Pc(V) 
of  armor/bullet  systems  was  always  monotonic. 

(2)  Direct  measurement  of  Pc(V)  in  AljOj/Kevlar  systems  confirmed  that  the  Pc(V)  was 
nonmonotonic. 

(3)  Observation  of  the  bullet-shattering  behavior  as  a  function  of  velocity  strongly  suggested 
that  the  nonmonotonic  Pc(V)  in  AljOj/Kevlar  systems  was  caused  by  the  shatter-gap 
behavior. 

(4)  A  proposed  V50  test  procedure  was  developed  applicable  for  testing  systems  with  either 
monotonic  or  nonmonotonic  Pc(V). 

(5)  Results  of  Phase  n  interlaboratory  V50  tests  program  using  the  proposed  Vjo  test  procedure 
showed  a  significantly  reduced  V50  data  scatter  as  compared  to  Phase  I  interlaboratory  data. 

(6)  To  further  reduce  other  sources  of  interlaboratory  Vjq  data  scatter,  we  recommend  frequent 
velocity  measurement  calibration  and  certification,  gun-noise  reduction,  bullet  pitch  and  yaw 
measurement  and  reduction,  and  standardization  of  the  target  fixture,  especially  for  ceramic/ 
composite  armor  plates. 
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Appendix  A: 

Proposed  V50  Test  Procedure 
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1 ,  Total  number  of  shots  -  15. 


2.  Target  fixture: 

For  cCTamic/composite  targets,  it  should  be  sandwiched  ti^tly  between  2  steel  plates  having  toe 
same  outside  dimensions  as  toe  target,  a  central  widow  at  least  2  in  larger  than  toe  ceramic  tile 
dimensions  and  a  contact  area  between  toe  steel  plates  and  toe  target  at  least  1.5  in  wide.  The 
target  should  be  mounted  to  a  rigid  holder  so  that  toe  bullet  impacts  at  0°  obliquity. 

3.  Use  a  0.020-in-toick  2024  aluminum  sheet,  ~6  in  behind  toe  armor,  as  a  witness  plate  to  record 
partial  or  complete  penetration. 

4.  Test  engineer  should  develop  and  use  his/her  own  powder  curve  for  bullets  to  control  projectile 
velocity  within  ±  25  fps. 

5.  Test  engineer  should  ensure  low  projectile  pitch-and-yaw  angle  (below  3  °)  for  toe  entire  velocity 
window  of  interest  and  use  yaw-card  or  flash  x-ray  to  record  pitch-and-yaw  for  each  shot. 
Discard  data  fi'om  shots  with  excessive  pitch-and-  yaw  (more  than  10°). 

6.  Test  OTgineer  should  determine  toe  velocity  window  of  interest:  between  Vjo*  to 

7.  For  low-  V50  test: 

number  of  shots  -  9,  Vj^.!  -  Vj  +  AV  (Vj  -  intended  velocities,  i  -  1-9) 

+  100  fps  (starting  velocity) 

AV  -  +100  fps  (if  previous  shot  is  partial) 

AV  -  - 100  ^s  (if  previous  shot  is  complete) 

8.  For  high- V50  test: 

number  of  shots  -  6,  Vj^.]  -  Vj  +  AV  (V;  -  intended  velocities,  i  -  10-15) 

V,o  -  Vygh  - 100  fps  (starting  velocity) 

AV  -  +100  fps  (if  previous  shot  is  partial) 

AV  -  - 100  ^s  (if  previous  shot  is  complete) 

9.  Use  a  statistical  computer  program,  MDMLE,  to  evaluate  toe  values  of  L0W-V50  and  High-Vjo 
firom  toe  firing  data. 

10.  Testing  engineer  should  prepare  firing  data  sheets  which  should  include  firing  sequence,  impact 
velocities,  partial/complete,  yaw,  Iow-Vjq  and  high-Vjo  obtained. 


Intentionally  left  blank. 
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Appendix  B: 

FORTRAN  Code  Listing  of  MDMLE  Program 
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c 

c  This  prograin  provides  a  Modified  Dixon  Maxunum  Likelihood  Estimate 

cforVSO. 

c 

program  mdmle 
character  flname*64 
dimension  stress(  100,2) 
c 

c  Establish  necessary  inputs, 
c 

callset(delta,nsample,stress4iumit,epsilonjionresp,flname) 

c 

c  Order  the  nonresponse  velocities  followed  by  the  response  velocities 
c 

call  sort(stress,nsample4ionresp) 
c 

c  Compute  the  Modified  Dixon  Method  MLE 
c 

call  mle(stress,v50,numit,epsilon4isample,delta4ionresp) 
c 

c  Display  results 
c 

call  out(delta,nsample, stress, v50) 
stop 
end 
c 

c  Subroutine  Set 
c 

c  This  subroutine  allows  for  data  to  be  read  in  from  a  file 
c  or  from  the  screen, 
c 

subroutine  set(delta4isample,stress4iumit,epsilon,nonresp,flname) 

dimension  stress(  100,2) 

character  flname*64,dansa*64,dansb*64 

nonresp-0 

write(*,*)  ‘Name  the  output  file’ 

read(*,’(a)’)flname 

write(*,*)  ‘Name  the  output  file’ 

read(*,’(a)’)flname 

open(8,file-flname,status=’new’) 

write(*,*)’Would  you  like  to  read  the  program  information?  y/n’ 

read(*,’(a)’)dansa 

if(dansa.eq.  ’y ’)then 

write(*,*)’This  FORTRAN  routine  delivers  a  maximum  likelihood’ 
write(*,*)’estimate  for  the  mean  of  the  logistic  response  distr-’ 
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write(*,*)’ibution.  This  is  a  suitable  substitute  for  the  normal’ 
write(*,*) ’distribution  over  the  .2  to  .8  quantiles.  Assumed  is’ 
write(*,*)’that  the  standard  deviation  of  the  response  distr-’ 
write(*,*)’ibution  is  known  to  be  the  “delta  velocity  change’’.’ 
write(*,*) 

write(*,*)’Because  of  this  assumption,  the  estimate  will  exist’ 
write(*,*)’if  at  least  one  nonresponse  and  one  response  are’ 
write(*,*)’present  in  the  data.  If  not,  a  code  of  -9999  will  be’ 
write(*,*) ’printed.  Occasionally,  the  Newton-Raphson  procedure’ 
write(*,*)’will  diverge.  A  code  of  -7777  is  given  when  this’ 
write(*,*) ’occurs.  If  convergence  has  not  been  achieved  in  the’ 
write(*,*)’number  of  iterations  allowed  by  the  user,  a  code  of 
write(*,*)’-8888  is  given.  Increase  the  number  of  iterations’ 
write(*,*)’allowed  or  relax  the  epsilon  for  convergence  to  ‘ 
write(*,*) ’correct  this.’ 
write(*,*) 

write(*,*)’Hit  enter  to  continue.’ 

read(*,’(a)’)dansb 

write(*,*) 

write(*,*)’There  is  no  attempt  on  the  part  of  the  author  to’ 
write(*,*)’pass  this  as  polished  software.  It  is  not.  It  is,’ 
write(*,*)’however,  functional  and  will  allow  the  user  to’ 
write(*,*) ’quickly  and  easily  get  a  maximum  likelihood  estimate’ 
write(*,*)’for  the  V50  instead  of  using  the  rough  approximation’ 
write(*,*)’that  involves  averaging  stimulus  levels.’ 
write(*,*) 

write(*,*)’Dr.  Barry  A.  Bodt’ 

write(*,*)’Army  Research  Lab’ 

write(*,*)’babodt@arl.army.mir 

write(*,*)’4  June,  1996’ 

writeC*,*) 

else 

endif 

write(*,*) 

write(*,*) ’Enter  the  delta  velocity  change.’ 
read(*,*)delta 

write(*,*)’Enter  the  number  of  number  of  shots  taken.’ 
read(*,*)nsample 

write(*,*) ’Enter  each  velocity,  followed  by  a  response  (1)’ 
write(*,*)’or  a  nonresponse  (0).’ 
do  10  j-l,nsample 
write(*,*)’test  ‘,j 
read(*,*)stress(j,  1  ),stress(j,2) 
if(stress(j,2).eq.0)nom-esp-nonresp+l 
10  continue 
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write(*,*) 

write(*,*)’How  many  iterations  are  to  be  allowed  for  the’ 
write(*,*)’Newton-Raphson  algorithm  and  what  is  the  epsilon’ 
write(*,*)’for  convergence?  Enter  the  number  of  iterations’ 
write(*,*)’followed  by  the  convergence  criterion.’ 
read(*,*)numit,epsilon 
return 
end 
c 

c  Subroutine  Sort 
c 

c  This  routine  is  just  a  simple  bubble  sort  to  order  the  data  in 
c  the  stress  array  according  to  increasing  velocity  within  response 
c 

subroutine  sort(stress4isample^onresp) 
dimmsion  stress(  100,2) 
c 

c  First  order  the  data  according  to  response  type  with 
c  nonresponse  (0)  and  response  (1). 
c 

do  30  i-l^isample-l 
do  20  j-i+1  ,nsample 
if(stress(j,2).lt.stress(i,2))then 
dolOk-1,2 
temp-stress(i,k) 
stress(idc)-stress(j,k) 
stress(j,k)-temp 
10  continue 
endif 

20  continue 
30  continue 
c 

c  Next  order  the  nonresponse  data  according  to  velocity, 
c 

do  60  ii-l,nonresp-l 
do  50ij-ii+l,nonresp 
if(stress(jj,  1 )  .lt.stress(ii,  1  ))then 
do40kk-l,2 
temp-stress(ii,kk) 
stress(ii,kk)-stress(jj3ck) 
stress(jj,kk)-temp 
40  continue 
endif 

50  continue 
60  continue 
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c  Finally  order  the  response  data  according  to  velocity, 
c 

do  90  iii-nonresp+1  ^sample-1 
do  80  jjj-iii+l  ^sample 
if(stress(iij,  1  ).lt.stress(iii,  1  ))then 
do70kkk-l,2 
temp-stress(iii4ddc) 
stress(iii4dck)-stress(jij[j4ckk) 
stress(jij4ckk)-temp 
70  continue 
endif 

80  continue 
90  continue 
return 
end 
c 

c  Output  of  initial  conditions 
c 

subroutine  out(delta4isample,stress,v50) 

dimension  stress(  100,2) 

write(8,*) 

write(*,*) 

write(8, 1 00)delta,nsample 

write(* ,  1 00)delta,nsample 

write(8,*) 

write(*,*) 

write(8,200) 

write(*,200) 

do  10  j-=l,nsample 

write(8,300)(stress(j,k),k»l  ,2) 

write(*,300)(stress(j,k),k-l  ,2) 

10  continue 
write(8,*) 
write(*,*) 
write(8,400)v50 
write(*,400)v50 
return 

100  formate  lx, ’delta  -  ‘,f8.3,5x,’shots  -  ‘,i8) 

200  format(lx,’velocity  response’) 

300  format(lx,f8.3,6x,f4.2) 

400  formate lx,’mle  v50  estimate  -  ‘,f8.3) 
end 

subroutine  mleestress,v50,numit,epsilon,nsample,delta,nonresp) 
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dimension  stress(  100,2) 
c 

c  check  to  see  that  the  estimate  exists  (at  least  one 
c  one  nonresponse  and  one  response 
c 

nresp-nsample-nonresp 

if(nonresp.eq.nsample.or.nresp.eq.nsample)then 

V50-9999 

return 

endif 

c 

c  With  the  logistic  model  f(alpha+beta*x),  the  standard 
c  deviation  is  related  to  beta  dirough  the  equation 
c  beta-pi/(sigma*sqrt(3)).  In  our  model,  the  standard 
c  deviation  is  being  taken  to  be  the  step  value  for  the 
c  data  collection,  hence  the  next  line  of  code, 
c 

beta-3.141592/(delta*sqrt(3)) 

c 

c  establish  a  reasonable  starting  value  for  the  Newton- 
c  Raphson  procedure  using  the  average  of  the  lowest  complete 
c  and  the  highest  partial.  Then  convert  that  v50  start  to 
c  the  parameters  for  the  f(alpha  +  beta*X)  formulation  to 
c  form  alphaO  for  the  Newton-Raphson. 
c 

xmuO-(stress(nonresp,  1  )+stress(nonresp+ 1 , 1  ))/2 
alphaO— beta*xmuO 
alphnml-alphaO 
do  20  k-l,numit 
xnumer-0 
xdenom-0 
c 

c  Compute  the  update  for  the  Newton  Raphson  approximation.  It 
c  involves  the  computation  of  probabilities  from  the  logistic 
c  distribution.  The  update  will  be  called  xmargin. 
c 

do  10  j-l,nsample 

dist- 1/(  1  +exp(-alphnm  1  -beta*  stress(j,  1 ))) 
xnumer-xnumer+stress(j  ,2)-dist 
xdenom-xdenom+dist*(  1  -dist) 

10  continue 
c 

c  Because  of  the  instability  of  the  Newton-Raphson  in  this 
c  situation,  it  is  possible  for  the  xdenom  above  to  be  very 
c  close  to  zero.  This  will  occur  when  the  “alpha  n  minus  1” 
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c  value,  the  current  parameter  estimate,  is  very  far  from 
c  the  actual  root.  Then  the  probabilites  will  be  computed 
c  for  velocity  levels  out  in  the  tail  of  the  distribution, 
c  Either  dist  or  1-dist  is  likely  to  produce  zeros  in  that 
c  instance.  A  good  starting  value  is  a  hedge  against  this, 
c  The  shape  of  the  function  for  which  the  root  is  desired 
c  is  such  Aat  if  you  are  out  a  little  too  far,  the  algorithm 
c  will  diverge,  not  converge.  This  is  likely  to  be  first 
c  manifested  in  a  zero  denominator, 
c 

if(xdenom.eq.O)then 
v50— 7777 
return 
endif 
c 

c  Update  the  parameter  alpha.  Check  to  see  if  it  satisfies  the 
c  convergence  criterion.  If  it  does,  compute  the  mean  based  on 
c  the  two  parameters  of  f(alpha  +  beta*X) 
c 

xmargin-xnumer/xdenom 

alphnml  -alphnml  +xmargin 

xtest-abs(xmargin) 

if(xtest.lt.epsilon)then 

alpha-alphnml 

v50— alpha/beta 

return 

endif 

20  continue 
c 

c  If  convergence  is  so  slow  that  it  cannot  occur  in  the  specified 
c  number  of  iterations,  assign  it  a  code  of  -8888.  This  is  unlikely 
c  to  occur,  but  is  possible, 
c 

v50— 8888 

return 

end 
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